SEV Met Yearly nlme Models

This report runs base, AR1, and AR2 models (variable ~ time) for annual SEV meteorological data using the nlme R package. Mean annual air temperature, minimum annual air temperature, maximum annual air temperature, and total annual precipitation are modeled. More sophisticated models will likely be required after this preliminary analysis.

---

library(tidyverse)
library(lubridate)
library(nlme)
library(MuMIn)
library(gridExtra)
library(kableExtra)

# Load yearly data ---------------------------------------------

path_to_files <- "../data/processed_data/"

# load and only keep vars of interest
y <- read_csv(paste0(path_to_files, "met_yearly_gap_filled.csv")) %>% 
  mutate(sta = as.factor(sta)) %>% 
  select(sta:ppt)

str(y)

## tibble [110 × 6] (S3: tbl_df/tbl/data.frame)
##  $ sta   : Factor w/ 4 levels "40","42","49",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ year  : num [1:110] 1988 1989 1990 1991 1992 ...
##  $ airt  : num [1:110] 13.8 13.8 13.8 12.5 12.3 ...
##  $ maxair: num [1:110] 35 35 35 37.7 38.6 ...
##  $ minair: num [1:110] -7.94 -7.94 -7.94 -13.76 -15.11 ...
##  $ ppt   : num [1:110] 151 140 232 335 241 ...

summary(y)

##  sta          year           airt           maxair          minair       
##  40:34   Min.   :1988   Min.   :11.91   Min.   :32.01   Min.   :-29.930  
##  42:33   1st Qu.:2001   1st Qu.:13.15   1st Qu.:37.13   1st Qu.:-16.495  
##  49:23   Median :2008   Median :13.82   Median :38.70   Median :-14.780  
##  50:20   Mean   :2007   Mean   :13.89   Mean   :38.50   Mean   :-14.900  
##          3rd Qu.:2015   3rd Qu.:14.67   3rd Qu.:39.84   3rd Qu.:-12.982  
##          Max.   :2021   Max.   :16.43   Max.   :48.72   Max.   : -3.797  
##       ppt       
##  Min.   :122.0  
##  1st Qu.:220.9  
##  Median :247.9  
##  Mean   :262.5  
##  3rd Qu.:301.7  
##  Max.   :599.4

---

Prepare data for modeling -

Add a time variable for nlme modeling and split data into individual met stations.

# function to select station and variable and add 'time' var
prepare_data_for_nlme <- function(data, sta) {
  # data = dataframe to use
  # sta = met station (in quotes)
  data %>% 
    filter(sta == {{ sta }}) %>% 
    mutate(time = as.numeric(as.factor(year)))
}


m40 <- prepare_data_for_nlme(y, "40")
m42 <- prepare_data_for_nlme(y, "42")
m49 <- prepare_data_for_nlme(y, "49")
m50 <- prepare_data_for_nlme(y, "50")

---

Set up functions -

This code produces several functions that will be useful for producing cleaner code.

# function for plotting preliminary yearly graph -----------------------------------------------------------

plot_yearly_prelim <- function(data, var, title, y_axis_label) {
  # data = dataset
  # var = variable (not in quotes)
  # title = title for graph (in quotes)
  # y_axis_label = label for y-axis (in quotes)
  ggplot(data, aes(x = year, y = {{ var }}, col = sta)) +
    geom_line(color = "burlywood") +
    geom_smooth(method = "lm", color = "black", size = 0.6) +
    facet_wrap(~ sta) +
    labs(title = title,
         x = "Year",
         y = y_axis_label)
}

# function for plotting graphs with best nlme model results ------------------------------------------------

plot_yearly_results <- function(data, var, coeffs, title, y_axis_label) {
  # data = dataset
  # var = variable (not in quotes)
  # coeffs = object containing coefficients and p-values of best model
  # title = title for graph (in quotes)
  # y_axis_label = label for y-axis (in quotes)
  ggplot(data, aes(x = year, y = {{ var }})) +
    geom_line(color = "burlywood") +
    geom_smooth(method = "lm", color = "black", size = 0.6) +
    labs(title = title,
         x = "Year",
         y = y_axis_label,
         caption = paste("slope = ", round(coeffs[2,1], 3), "\n(p-value = ", coeffs[2,2], ")"))
}


# functions to run nlme models ----------------------------------------------------------

# base model function
base_nlme_model <- function(data, var) {
  # data = dataframe
  # var = var to model
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      method = "ML",
      na.action = na.omit)
}



# AR1 model function
AR1_nlme_model <- function(data, var) {
  # data = dataframe
  # var = var to model
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      correlation = corARMA(form = ~ time, p = 1),
      method = "ML",
      na.action = na.omit)
}


# AR2 model function
AR2_nlme_model <- function(data, var) {
  # data = dataframe
  # var = var to model
  
  model_specification = as.formula(paste0(var, " ~ time"))
  
  gls(model = model_specification,
      data = data,
      correlation = corARMA(form = ~ time, p = 2),
      method = "ML",
      na.action = na.omit)
}

---

Mean Annual Air Temperature (C) -

plot_yearly_prelim(y, airt, "Annual Mean Air Temperature", "Temperature (C)")

## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "airt")
m40_mAR1 <- AR1_nlme_model(m40, "airt")
m40_mAR2 <- AR2_nlme_model(m40, "airt")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m40_mbase	3	53.68867
m40_mAR1	4	55.62003
m40_mAR2	5	58.37968

summary(m40_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   52.88867 57.46775 -23.44433
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 13.271801 0.17431449 76.13711  0.0000
## time         0.021443 0.00868861  2.46794  0.0191
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1560729 -0.5446061  0.2023188  0.6881964  1.6138251 
## 
## Residual standard error: 0.4821986 
## Degrees of freedom: 34 total; 32 residual

m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]

(m40_plot <- plot_yearly_results(m40, airt, m40_coeffs, "Met 40 - Deep Well \nAnnual Mean Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "airt")
m42_mAR1 <- AR1_nlme_model(m42, "airt")
m42_mAR2 <- AR2_nlme_model(m42, "airt")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m42_mbase	3	50.38362
m42_mAR1	4	52.94158
m42_mAR2	5	55.65049

summary(m42_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   49.55604 54.04556 -21.77802
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 12.618607 0.17205515 73.34048  0.0000
## time         0.011535 0.00883012  1.30629  0.2011
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -1.987976067 -0.644320891 -0.002985058  0.916472339  1.576226874 
## 
## Residual standard error: 0.468135 
## Degrees of freedom: 33 total; 31 residual

m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]

(m42_plot <- plot_yearly_results(m42, airt, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Mean Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "airt")
m49_mAR1 <- AR1_nlme_model(m49, "airt")
m49_mAR2 <- AR2_nlme_model(m49, "airt")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m49_mbase	3	34.60607
m49_mAR1	4	37.16381
m49_mAR2	5	39.81788

summary(m49_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   33.34291 36.74939 -13.67145
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 14.203751 0.19776763 71.82041  0.0000
## time         0.024676 0.01442369  1.71081  0.1018
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1832629 -0.6706411  0.2109746  0.6989873  1.7965345 
## 
## Residual standard error: 0.4384421 
## Degrees of freedom: 23 total; 21 residual

m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]

(m49_plot <- plot_yearly_results(m49, airt, m49_coeffs, "Met 49 - Five Points \nAnnual Mean Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "airt")
m50_mAR1 <- AR1_nlme_model(m50, "airt")
m50_mAR2 <- AR2_nlme_model(m50, "airt")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m50_mbase	3	33.37573
m50_mAR1	4	35.87783
m50_mAR2	5	39.37225

summary(m50_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   31.87573 34.86293 -12.93786
## 
## Coefficients:
##                 Value  Std.Error  t-value p-value
## (Intercept) 14.831383 0.22625535 65.55152  0.0000
## time         0.051514 0.01888743  2.72741  0.0138
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.81172625 -0.64841217 -0.09216346  0.63402578  1.96854695 
## 
## Residual standard error: 0.462067 
## Degrees of freedom: 20 total; 18 residual

m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]

(m50_plot <- plot_yearly_results(m50, airt, m50_coeffs, "Met 50 - Blue Grama \nAnnual Mean Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Mean Annual Air Temperature

# airt results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

---

Minimum Annual Air Temperature (C) -

plot_yearly_prelim(y, minair, "Annual Minimum Air Temperature", "Temperature (C)")

## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "minair")
m40_mAR1 <- AR1_nlme_model(m40, "minair")
m40_mAR2 <- AR2_nlme_model(m40, "minair")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m40_mbase	3	189.6589
m40_mAR1	4	191.1119
m40_mAR2	5	193.7564

summary(m40_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC     BIC    logLik
##   188.8589 193.438 -91.42944
## 
## Coefficients:
##                  Value Std.Error    t-value p-value
## (Intercept) -13.656159 1.2874555 -10.607092  0.0000
## time         -0.136817 0.0641725  -2.132025  0.0408
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.6474716 -0.3773667  0.0234212  0.4512047  1.7215955 
## 
## Residual standard error: 3.561432 
## Degrees of freedom: 34 total; 32 residual

m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]

(m40_plot <- plot_yearly_results(m40, minair, m40_coeffs, "Met 40 - Deep Well \nAnnual Minimum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "minair")
m42_mAR1 <- AR1_nlme_model(m42, "minair")
m42_mAR2 <- AR2_nlme_model(m42, "minair")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m42_mbase	3	189.4901
m42_mAR1	4	190.6799
m42_mAR2	5	192.9317

summary(m42_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC     BIC    logLik
##   188.6625 193.152 -91.33124
## 
## Coefficients:
##                  Value Std.Error   t-value p-value
## (Intercept) -13.220095 1.4158546 -9.337184  0.0000
## time         -0.094695 0.0726637 -1.303196  0.2021
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -3.52824489 -0.32599081 -0.05008804  0.56078847  2.47065756 
## 
## Residual standard error: 3.852318 
## Degrees of freedom: 33 total; 31 residual

m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]

(m42_plot <- plot_yearly_results(m42, minair, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Minimum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "minair")
m49_mAR1 <- AR1_nlme_model(m49, "minair")
m49_mAR2 <- AR2_nlme_model(m49, "minair")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m49_mbase	3	123.2896
m49_mAR1	4	125.6878
m49_mAR2	5	128.9351

summary(m49_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   122.0264 125.4329 -58.01321
## 
## Coefficients:
##                  Value Std.Error    t-value p-value
## (Intercept) -14.590316 1.3596668 -10.730802    0.00
## time          0.006294 0.0991639   0.063475    0.95
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -4.0179895 -0.3267865  0.1783908  0.5709077  1.3371903 
## 
## Residual standard error: 3.014321 
## Degrees of freedom: 23 total; 21 residual

m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]

(m49_plot <- plot_yearly_results(m49, minair, m49_coeffs, "Met 49 - Five Points \nAnnual Minimum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "minair")
m50_mAR1 <- AR1_nlme_model(m50, "minair")
m50_mAR2 <- AR2_nlme_model(m50, "minair")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m50_mbase	3	112.1543
m50_mAR1	4	115.2297
m50_mAR2	5	118.1863

summary(m50_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   110.6543 113.6415 -52.32717
## 
## Coefficients:
##                  Value Std.Error   t-value p-value
## (Intercept) -14.264158 1.6215367 -8.796691  0.0000
## time          0.072444 0.1353633  0.535179  0.5991
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8049348 -0.3200412  0.1852600  0.5306493  1.1520424 
## 
## Residual standard error: 3.311562 
## Degrees of freedom: 20 total; 18 residual

m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]

(m50_plot <- plot_yearly_results(m50, minair, m50_coeffs, "Met 50 - Blue Grama \nAnnual Minimum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Minimum Annual Air Temperature

# minair results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

---

Maximum Annual Air Temperature (C) -

plot_yearly_prelim(y, maxair, "Annual Maximum Air Temperature", "Temperature (C)")

## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "maxair")
m40_mAR1 <- AR1_nlme_model(m40, "maxair")
m40_mAR2 <- AR2_nlme_model(m40, "maxair")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m40_mbase	3	140.3235
m40_mAR1	4	135.6487
m40_mAR2	5	136.8874

The AR1 model has the lowest AICc score. This model is used for plotting.

summary(m40_mAR1)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   134.2694 140.3748 -63.13468
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##       Phi 
## 0.4530549 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 37.64031 0.9617813 39.13604  0.0000
## time         0.08334 0.0473790  1.75903  0.0881
## 
##  Correlation: 
##      (Intr)
## time -0.862
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6634576 -0.5746671 -0.1493423  0.6655313  3.1746966 
## 
## Residual standard error: 1.732308 
## Degrees of freedom: 34 total; 32 residual

m40_coeffs <- summary(m40_mAR1)$tTable[,c(1, 4)]

(m40_plot <- plot_yearly_results(m40, maxair, m40_coeffs, "Met 40 - Deep Well \nAnnual Maximum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "maxair")
m42_mAR1 <- AR1_nlme_model(m42, "maxair")
m42_mAR2 <- AR2_nlme_model(m42, "maxair")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m42_mbase	3	129.8644
m42_mAR1	4	127.7573
m42_mAR2	5	130.0821

The AR1 model has the lowest AICc score. This model is used for plotting.

summary(m42_mAR1)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   126.3287 132.3147 -59.16435
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##       Phi 
## 0.3867261 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 34.94019 0.8290183 42.14647  0.0000
## time         0.06777 0.0421570  1.60763  0.1181
## 
##  Correlation: 
##      (Intr)
## time -0.864
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.94732099 -0.54279907 -0.08489296  0.66776673  2.67346131 
## 
## Residual standard error: 1.572183 
## Degrees of freedom: 33 total; 31 residual

m42_coeffs <- summary(m42_mAR1)$tTable[,c(1, 4)]

(m42_plot <- plot_yearly_results(m42, maxair, m42_coeffs, "Met 42 - Cerro Montoso \nAnnual Maximum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "maxair")
m49_mAR1 <- AR1_nlme_model(m49, "maxair")
m49_mAR2 <- AR2_nlme_model(m49, "maxair")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m49_mbase	3	116.1389
m49_mAR1	4	117.1260
m49_mAR2	5	116.5718

The base model has the lowest AICc and is used for plotting.

summary(m49_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   114.8758 118.2822 -54.43788
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 41.07945 1.1639166 35.29415  0.0000
## time        -0.10151 0.0848873 -1.19584  0.2451
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.3311794 -0.5068364 -0.1698531  0.1037714  3.2757736 
## 
## Residual standard error: 2.580352 
## Degrees of freedom: 23 total; 21 residual

m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]

(m49_plot <- plot_yearly_results(m49, maxair, m49_coeffs, "Met 49 - Five Points \nAnnual Maximum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "maxair")
m50_mAR1 <- AR1_nlme_model(m50, "maxair")
m50_mAR2 <- AR2_nlme_model(m50, "maxair")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m50_mbase	3	63.96627
m50_mAR1	4	63.49305
m50_mAR2	5	66.10862

The AR1 model has the lowest AICc score. This model is used for plotting.

summary(m50_mAR1)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   60.82638 64.80931 -26.41319
## 
## Correlation Structure: AR(1)
##  Formula: ~time 
##  Parameter estimate(s):
##        Phi 
## -0.4018488 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 38.90612 0.3271151 118.93708  0.0000
## time         0.07382 0.0274879   2.68568  0.0151
## 
##  Correlation: 
##      (Intr)
## time -0.882
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.65156196 -0.77367285 -0.05306683  0.49518631  2.02742794 
## 
## Residual standard error: 0.9854843 
## Degrees of freedom: 20 total; 18 residual

m50_coeffs <- summary(m50_mAR1)$tTable[,c(1, 4)]

(m50_plot <- plot_yearly_results(m50, maxair, m50_coeffs, "Met 50 - Blue Grama \nAnnual Maximum Air Temperature", "Temperature (C)"))

## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Maximum Annual Air Temperature

# maxair results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

---

Total Annual Precipitation (mm) -

plot_yearly_prelim(y, ppt, "Total Annual Precipitation", "Precipitation (mm)")

## `geom_smooth()` using formula 'y ~ x'

Met 40:

# Met 40:

m40_mbase <- base_nlme_model(m40, "ppt")
m40_mAR1 <- AR1_nlme_model(m40, "ppt")
m40_mAR2 <- AR2_nlme_model(m40, "ppt")

AICc scores:

AICc(m40_mbase, m40_mAR1, m40_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m40_mbase	3	376.9800
m40_mAR1	4	379.0143
m40_mAR2	5	380.8387

summary(m40_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##      AIC      BIC  logLik
##   376.18 380.7591 -185.09
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 239.73384 20.234522 11.847764  0.0000
## time         -0.79358  1.008579 -0.786827  0.4372
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.7633224 -0.6843958  0.0941523  0.8117423  1.7513612 
## 
## Residual standard error: 55.97387 
## Degrees of freedom: 34 total; 32 residual

m40_coeffs <- summary(m40_mbase)$tTable[,c(1, 4)]

(m40_plot <- plot_yearly_results(m40, ppt, m40_coeffs, "Met 40 - Deep Well \nTotal Annual Precipitation", "Precipitation (mm)"))

## `geom_smooth()` using formula 'y ~ x'

Met 42:

# Met 42:

m42_mbase <- base_nlme_model(m42, "ppt")
m42_mAR1 <- AR1_nlme_model(m42, "ppt")
m42_mAR2 <- AR2_nlme_model(m42, "ppt")

AICc scores:

AICc(m42_mbase, m42_mAR1, m42_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m42_mbase	3	386.6521
m42_mAR1	4	389.0729
m42_mAR2	5	390.7934

summary(m42_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC     BIC    logLik
##   385.8245 390.314 -189.9122
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 391.6703 28.079405 13.948670  0.0000
## time         -3.5405  1.441075 -2.456826  0.0198
## 
##  Correlation: 
##      (Intr)
## time -0.872
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.60096602 -0.53817435 -0.06228751  0.56988895  3.13557611 
## 
## Residual standard error: 76.39964 
## Degrees of freedom: 33 total; 31 residual

m42_coeffs <- summary(m42_mbase)$tTable[,c(1, 4)]

(m42_plot <- plot_yearly_results(m42, ppt, m42_coeffs, "Met 42 - Cerro Montoso \nTotal Annual Precipitation", "Precipitation (mm)"))

## `geom_smooth()` using formula 'y ~ x'

Met 49:

# Met 49:

m49_mbase <- base_nlme_model(m49, "ppt")
m49_mAR1 <- AR1_nlme_model(m49, "ppt")
m49_mAR2 <- AR2_nlme_model(m49, "ppt")

AICc scores:

AICc(m49_mbase, m49_mAR1, m49_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m49_mbase	3	253.8175
m49_mAR1	4	256.5620
m49_mAR2	5	259.8661

summary(m49_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC    logLik
##   252.5543 255.9608 -123.2772
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 254.18442 23.215104 10.949096  0.0000
## time         -1.99381  1.693136 -1.177584  0.2521
## 
##  Correlation: 
##      (Intr)
## time -0.875
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.01258367 -0.67241549  0.07653863  0.62248810  2.09072501 
## 
## Residual standard error: 51.46686 
## Degrees of freedom: 23 total; 21 residual

m49_coeffs <- summary(m49_mbase)$tTable[,c(1, 4)]

(m49_plot <- plot_yearly_results(m49, ppt, m49_coeffs, "Met 49 - Five Points \nTotal Annual Precipitation", "Precipitation (mm)"))

## `geom_smooth()` using formula 'y ~ x'

Met 50:

# Met 50:

m50_mbase <- base_nlme_model(m50, "ppt")
m50_mAR1 <- AR1_nlme_model(m50, "ppt")
m50_mAR2 <- AR2_nlme_model(m50, "ppt")

AICc scores:

AICc(m50_mbase, m50_mAR1, m50_mAR2) %>% 
  kbl() %>% 
  kable_styling(full_width = FALSE, position = "left")

	df	AICc
m50_mbase	3	222.6721
m50_mAR1	4	225.6402
m50_mAR2	5	229.2037

summary(m50_mbase)

## Generalized least squares fit by maximum likelihood
##   Model: model_specification 
##   Data: data 
##        AIC      BIC   logLik
##   221.1721 224.1593 -107.586
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 245.26421 25.695541 9.545011  0.0000
## time          0.26436  2.145022 0.123244  0.9033
## 
##  Correlation: 
##      (Intr)
## time -0.877
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.95082391 -0.49761448 -0.02635629  0.88260113  1.46242374 
## 
## Residual standard error: 52.47638 
## Degrees of freedom: 20 total; 18 residual

m50_coeffs <- summary(m50_mbase)$tTable[,c(1, 4)]

(m50_plot <- plot_yearly_results(m50, ppt, m50_coeffs, "Met 50 - Blue Grama \nTotal Annual Precipitation", "Precipitation (mm)"))

## `geom_smooth()` using formula 'y ~ x'

Summary Plots of Total Annual Precipitation

# ppt results 
grid.arrange(m40_plot, m42_plot, m49_plot, m50_plot, ncol=2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'